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It is shown that the heart rate variability (HRV) in healthy and diseased humans possesses a 
hierarchical structure of the She-Leveque (SL) form. This structure, first found in measurements 
in turbulent fluid flows, implies further details in the HRV multifractal scaling. The potential of 
diagnosis is also discussed based on the characteristics derived from the SL hierarchy. 

PACS numbers: 87.19.Hh, 87.10.+e 

The heart beat interval in humans is known to exhibit fluctuation which is referred to as heart rate variability 
(HRV). Power spectrum analysis of the fluctuation revealed a 1 / /-like scaling pj. Recent studies indicated that 
healthy human HRV exhibits even higher complexity which can be characterized by multifractal scaling 0, Q ■ I n 
contrast, HRV in the pathological state such as congestive heart failure exhibits more monofractal-like scaling 0. 
The change of the HRV 1/f law in congestive heart failure is consistent with this result 4]. Such a multifractal- 
monofractal transition was also reported in parasympathetic nervous system (PNS) blockade experiment Q. Hence 
the manifestation of multifractal HRV is indicative of the proper autonomic regulation of the heart rate. Further 
studies revealed that the multifractal HRV have properties analogous to those found in fluid turbulence \6]. However, 
there is little understanding beyond the phenomenological description of multifractal HRV. 

In this paper, we exploit further the analogy of HRV to fluid turbulence and show the existence of a hierarchical 
structure in healthy and diseased HRV. This structure allows us to model the multifractality of HRV and make 
conjecture to the heart beat dynamics responsible for the multifractal scaling. The hierarchy, first proposed by She 
and Leveque (SL) to understand the statistical properties of turbulent flows, provides a successful framework to discuss 
and characterize the deviation from Kolmogorov monofractal scaling in fluid turbulence ||- When applied to study 
HRV, the SL hierarchy provides a model structure which possesses two advantages: (a) it simplifies the functional 
description of the multiscaling by using a maximum of only three parameters, and (b) it contains predictive power 
for HRV scaling in physiological states related to PNS withdrawal. One immediate implication is the potential use of 
this notion in applications such as diagnosis. 

Let the beat-to-beat RR interval (RRi) be r(i), where t is the beat number, and its increment be Ar(r) = r(t + 
t) — r(i). The SL hierarchy implies, for a range of r, 
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Here < f3 < 1 is a parameter of the hierarchy, A p , a function of p, S p (t) = (|Ar(r)| p ), the p-th order moment of 
|Ar(r)| denoted as the p-th order RRi structure function, 5°°(t) = linip^oo S p+ i(t) / S p (t) and (•) denotes statistical 
average. Since S°°(r) is dominated by the statistics of large Ar(r), it characterizes the most intense fluctuation in 
HRV. Moreover, given the empirical law S p (t) ~ in HRV the hierarchy implies |(| the scaling model 
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where ho and C are two other parameters of the hierarchy. It follows from (|TJ and J2J that 5°° (r) ~ t h ° . A nonlinear 
functional dependence of on p indicates multifractal scaling. Thus, the parameter (3 measures the degree of 
multifractality. In particular, (3 — > 1 leads to monofractal scaling. In the multifractal description of fluid turbulence, 
the parameter C can be shown to be the codimension of the most intense structure of the flow Since it is not yet 
possible to write down the equation of motion for long-term cardiovascular dynamical system, we assume a working 
definition for C as the "codimension parameter" of the signal. 

We follow the procedure developed in Ref. Q to check whether the RRi data possess a SL hierarchical form. This 
approach, based on the scaling property implied by the hierarchy, the so-called generalized extended self-similarity 
(GESS) in fluid turbulence j describes a power-law relationship between the normalized structure functions: 
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In the case of SL hierarchy, the exponents p n (p, q) depends only on the model parameter (3: 
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It follows that 

A Pn{p+ 5 P , q) = ^A PnM - ^^-f^ (5) 

where Ap n (p, q) = p n (p + dp, q) — p n (p, <?)• One can then plot Ap n (p + Sp, q) vs Ap n (p, q) to check (|5jl and hence the 
validity of the SL hierarchy. We use several databases to perform the calculations (3) and (5). The first database 
(DB1) contains 10 sets of daytime ambulatory RRi recordings taken from healthy young adults 5]. The second 
database (DB2) contains 18 sets of daytime normal sinus rhythm RRi data downloaded from public domain ^(J- We 
also analyze RRi data from congestive heart failure patients (DBCHF) from the same public domain [Tjj to study 
the intrincity of the hierarchy. 

Except for a few cases where excessive ectopic beats in the RRi data complicates the calculation of S p (t), and were 
therefore discarded from the analysis, GESS is found in both healthy and CHF HRV. The exponent p n {p,q) is then 
estimated from © and used in 10 to calculate Ap n (p,q). Typical Ap n {p + Sp,q) vs Ap n (p,q) plots are shown in 
Fig. 1. The observed linear trend implies JSJ. Hence, SL hierarchy is compatible with the multifractal scaling in HRV. 
From such plots, we estimate the value of j3 simultaneously from the slope and the intercept of the fitted straight 
lines. The results are given in Fig. 2a. The /3's from healthy HRV (DB1,DB2) cluster in the range [0.65,0.85] with 
those C(p) showing less curvature being characterized by larger (3 values (Fig. 1). The ft's from DBCHF are generally 
larger in values due to the monofractal-like scaling. 

To gain insight of the hierarchy, She and Waymire (SW) arrived at the hierarchy Q using multiplicative random 
cascade [Tl| . Their cascade consists of two dynamic components. One is the basic component that generates the 
singular dynamics over a continuum of scales. It can be shown that this dynamical component gives rise to the scaling 
term hop in J2J. SW's cascade contains an extra component, which they called the "defect dynamics" (DD), that 
modulates the singular structure through the multiplication of [3 in discrete steps |Tl) . It can be shown that this 
dynamical component contributes to the nonlinear term C(l — (3 P ) in J5J. To apply SW's interpretation requires some 
clarification since continuous scale invariance does not exist in HRV, due to the fact that RRi fluctuation between 
heart beats cannot be defined. This suggests a dominant DD in the generation of the multifractal scaling of HRV and 
a scaling model with ho ~ 0. Equivalently, this implies a hierarchy with a r-independent S°°. Since S°° cannot be 
directly calculated, to verify such a model we first re-write JTJ as 
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where p(p,q) = (1 — (3 P )/{1 — (3 q ) 0. We then form the quotient of iJBJ at distinct values of r and To, which, after 
some algebra, yields 
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Hence, F p ^(t, to) is independent of p and q and a r-independent S°°(t) implies a "constant" F p ^ q (r, tq) over a range 
of t and To values. Figure 3 shows F Ptq (T,To) for healthy and CHF HRV. It is seen that the condition ho ~ can 
be statistically ascertained and that ((p) ~ C(l — (3 P ). Given this, C is obtained by averaging £(p)/(l — (3 P ) over 
a range of p and its result has been shown in Fig. 2b. Moreover, C, as a function j3, shows an increasing trend as 
j3 — * 1 (Fig. 4). This functional relationship is consistent with the observations ho ~ and that £(p) becomes almost 
proportional to p as (3 —* 1. It can also be inferred from the earlier experimental studies, as we now explain. 

Recall that the fractal dimension D(h) of the set with a local scaling exponent h is related to ((p) through a 
Legendre transform: 

D(h) = min p \ph + d-((p)] (8) 
where d is the dimension of the embedding space. From J5| and with ho ~ 0, D(h) can be explicitly obtained as: 
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Let h* be the scaling exponent of the singular set with largest dimension, i.e., D(h*) is the maximum. For HRV, h* 
was found to increase its value from the multifractal-like scaling in healthy state to the monofractal-like scaling in the 
diseased and pathological states 0, Q ■ Using © , h* is derived explicitly as 



h* = C\n(l/(3). 
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If C is constant, h* decreases as monofractal scaling is approached (f3 — > 1), which contradicts what was observed 
In order for h* to increase as (3 — > 1, C{(3) must diverge as 1/ ln(l//3). Indeed, we find that the dependence of C 
on /3 can be well described by 0.2/3/ ln(l//3) (Fig. 4). Thus we have the result h* ~ 0.2/3 with ft,* increasing with (3 in 
accord with the experimental observations 0L| Finally, the asymptotic behaviour of C (as — » 1) presents a more 
favorable condition for diagnosis: i.e., the more monofractal- like scaling falls into the range with larger "separation" 
of the C values (see also Fig. 2b). 

In summary, we show that a hierarchical structure of the SL form exists in the healthy and diseased HRV. This 
property allows us to model the multifractal HRV in terms of only two parameters C and /3. Interestingly, C and /3 
are related by an empirical law captured in Fig. 4. This finding is important for two reasons. First, the empirical law 
appears universal and is capable of describing both healthy and CHF data. Second, the divergence of C as /3 — > 1 
implies potential in diagnosis using the hierarchical structure. To find a model that is compatiable with the current 
finding, we adopted SW's cascade which leads to further implications beyond the phcnomenological description of 
multifractal HRV scaling. This model is appealing in that it contains a modulating component (DD) acting on the 
singular dynamics to "tame" the fluctuation (since (3 < 1.) Its effect reminds us of the function of feedback regulation 
in biological systems. Hence, the SW's model provides a concrete example of how (additive) feedback mechanism 
may be integrated multiplicatively in a cascading structure to produce the observed HRV phenomenology. Further 
experiments will be needed to test and quantify these possibilities in more detailed physiological terms. 
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Fig. 1 

FIG. 1: Typical Ap n (p + Sp,q) vs. Ap n (p,q) plots for (a) healthy HRV and (b) CHF HRV with Sp = 0.2. We used q = 1 
and n = 2 for healthy HRV and q — 0.8 and n = 1.2 for CHF HRV. The arrow direction indicates the increasing p direction. 
The estimated slopes for the shown cases are f3 Sp = 0.9422 and 0.9881, respectively. Hence, j3 ~ (0.9422) 5 = 0.7425 and 
(0.9881) 5 ~ 0.9491, respectively. The insets show the corresponding £(p). It is seen that a large [3 in CHF HRV implies a ((p) 
with less curvature, i.e., monofractal-like scaling. 
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FIG. 2: Estimated values for model parameters (a) /3 and (b) C. Estimates from DB1, DB2 for 24 healthy subjects are given 
in circles and those from DBCHF for 30 congestive heart failure subjects are given in squares. Means and standard deviations 
in both cases are also shown by the bar-chart. 
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FIG. 4: C vs. j3 empirical law for healthy subjects (open circles) and congestive heart failure patients (solid circles). The solid 
line is the fit C = 0.2/3/ ln(l//3). 



